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Abstract. The present work considers diffusive shock ac- 
celeration at non-relativistic shocks using a system of 
stochastic differential equations (SDE) equivalent to the 
Fokker-Planck equation. We compute approximate solu- 
tions of the transport of cosmic particles at shock fronts 
with a SDE numerical scheme. Only the first order Fermi 
process is considered. The momentum gain is given by 
implicit calculations of the fluid velocity gradients using 
a linear interpolation between two consecutive time steps. 
We first validate our procedure in the case of single shock 
acceleration and retrieve previous analytical derivations of 
the spectral index for different values of the Peclet num- 
ber. The spectral steepening effect by synchrotron losses is 
also reproduced. A comparative discussion of implicit and 
explicit schemes for different shock thickness shows that 
implicit calculations extend the range of applicability of 
SDE schemes to infinitely thin ID shocks. The method is 
then applied to multiple shock acceleration that can be 
relevant for Blazar jets and accretion disks and for galac- 
tic centre sources. We only consider a system of identical 
shocks which free parameters are the distance between 
two consecutive shocks, the synchrotron losses time and 
the escape time of the particles. The stationary distribu- 
tion reproduces quite well the flat differential logarithm 
energy distribution produced by multiple shock effect, and 
also the piling-up effect due synchrotron losses at a mo- 
mentum where they equilibrate the acceleration rate. At 
higher momenta particle losses dominate and the spec- 
trum drops. The competition between acceleration and 
loss effects leads to a pile-up shaped distribution which ap- 
pears to be effective only in a restrict range of inter-shock 
distances of ~ 10-100 diffusion lengths. We finally com- 
pute the optically thin synchrotron spectrum produced 
such periodic pattern which can explain flat and/or in- 
verted spectra observed in Flat Radio spectrum Quasars 
and in the galactic centre. 



Key words: acceleration of particles - shock waves - radi- 
ation mechanisms - galaxy: centre - BL Lacertae objects 



1. Introduction 

The non-thermal radio and high-energy radiation spec- 
tra from extragalactic and galactic objects such as X-ray 
binaries, micro-Quasars, active galactic nuclei, jets and 
gamma-ray bursts require the presence of non-thermal 
particle distributions. One of the prime mechanisms for 
producing such distributions is diffusive acceleration at 
a shock front. This theory assumes the particle distribu- 
tions are isotropised by efficient scattering on wave tur- 
bulence on both sides of the shock and gain energy by 
a first order Fermi process upon crossing it. The equa- 
tion governing the particle distribution is of the Fokker- 
Planck type, containing both advection (dynamical fric- 
tion) and diffusion terms in the spatial variables and an 
advection term in the energy variable (or in the magni- 
tude of the particle momentum p), which is proportional 
to the divergence of the fluid flow, and remains valid in 
an integral sense even across shock fronts. Many gener- 
alisations of this equation have been discussed, but for 
the application to non-thermal radiation spectra, the most 
relevant are the inclusion of synchrotron losses (Webb et 
al 1984) and the extension to systems containing multi- 



ple shock fronts (Blandford fc Ost riker [1980] Spr uit [1988 , 
Achterber g |1990| , Schneider pJ93| Melrose |199l| Melrose 
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Crouch 1997] ). However, because of the difficulty of solv- 
ing the Fokker-Planck equation when the coefficients are 
complicated (and possibly discontinuous) functions of po- 
sition, these papers either adopted an idealised situation, 
or developed approximation schemes valid in only part of 
the parameter space. A particularly interesting alternative 
approach is to use the equivalence of the Fokker-Planck 
equation to a system of stochastic differential equations 
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(SDE's). Numerical integration of these is then akin to a 
Monte-Carlo simulation of the problem, which is relatively 
simple to implement, applies to complex flow patterns and 
places no restriction on the number of phase-space dimen- 
sions (e.g., Gardiner 1983). Several conventional Monte- 



Carlo simulations of part icle a cceleration exist (for a re- 
view see Jones & Ellison 1991), usually assuming a pre- 
scribed form of the mean free path of the particle as a 
function of the particle rigidity and plasma density. They 
have the advantage of being able to describe the evolution 
of both the thermal and suprathcrmal populations, as well 
as the non-linear back reaction of the non-thermal compo- 
nent on the shock profile. Relativistic shocks and large an- 
gle scattering are also e asy to include in such simulations 
(see Ellison et al. 1990). However, this approach always 



makes additional assumptions concerning the particle tra- 
jectory (e.g., that it is unperturbed between isotropising 
'collisions') which go beyond the diffusion approximation. 
In comparison, the SDE approach adopts the diffusion ap- 
proximation for test particle and is restricted to the trans- 
port of suprathermal particles. 

It is only recently that SDE systems have been ap- 
plied to astrophysical problems: in solar physics with the 
investigation of acceleration of fast electrons in the solar 

Conway et al 



corona (MacKinnon & Craig 1991 



1998) 



and in space physics with ion acceleration at the solar ter- 
mination shock (Chalov et al. 1995| ). In two papers Achter 
berg & Kriills Ql992| ) and Kriills & Achterberg ( [1994 - 



henceforth KA94), have applied the SDE approach to the 
problem of particle acceleration at astrophysical shocks, 
including the possibility of second-order Fermi accelera- 
tion (Schlickeiser 1989 ). These authors drew several im- 
portant conclusions. In particular, they showed that for 
Kolmogorov turbulence the second-order acceleration ef- 
fect is restricted to a small momentum range close to 
that of injection, because of the increase of the diffusion 
time with momentum. At high momentum, in contrast, 
the spectrum is formed by the competing effects of syn- 
chrotron losses and the first-order Fermi process. 

The numerical scheme used in KA94 is explicit, that 
is, the first-order Fermi acceleration term in the SDE is 
integrated forwards in time using the value of the fluid ve- 
locity gradient at the beginning of each step. This method 
requires a step short enough to resolve sharp features in 
the flow, such as a shock transition, and thus severely 
limits the ability to simulate acceleration in complex flow 
patterns containing structure on a large range of spatial 
scales. Our purpose in this paper is to propose and test an 
implicit method of integrating the SDE's. By implicit we 
mean that the coefficient of the first-order Fermi accelera- 
tion term is computed by linear interpolation between the 
end points of each time step. Such a scheme allows one 
to treat discontinuous shock structures using a finite time 
step and thus opens up the prospect of finding approxi- 
mate numerical solutions of the advection-diffusion equa- 
tion in complex flow patterns where the gradients of the 



fluid flow may become large over a distance much shorter 
than the shock thickness. This is our longer term goal; 
for the present, we limit ourselves to the problem of a ID 
single shock or to systems of multiple shocks. 

The organisation of the paper is as follows: in Sect. |^ 
we present the SDE system equivalent to the advection- 
diffusion equation and describe the explicit and implicit 
integration methods. We then test the implicit scheme in 
Sect. H by applying it to the (well-known) problem of ac- 
celeration at a single isolated shock front, both with and 
without synchrotron losses. We also compare the perfor- 
mances of implicit and explicit schemes. Section || focuses 
on the application to acceleration in a system of multiple 
shocks. Here, as well as testing the code against known 
approximate solutions, we present new solutions valid in 
regions of parameter space inaccessible to the analytic 
methods. We use the results to calculate the optically thin 
synchrotron radiation produced by the stationary particle 
distribution within a periodic pattern of shocks. In Sect. || 
we summarise our results and discuss extensions of the 
method to more complex problems. 

2. Monte-Carlo simulations 

2.1. Formulation of the SDE's 

The usual form of the advection-diffusion e quatio n, de- 
scribing the transport of cosmic-rays (Skilling 1975 ) is (in 
3D) 



f + (u.V)/-i(V.«)pg = S. 



D. 



(1) 



where S is a source term and D is a diffusion operator of 
the form 



D = V (TsV/) 



(2) 



The diffusion tensor 7t describes the spatial transport of 
particles. Adding the effects of synchrotron losses, and 
second-order Fermi acceleration one finds 



D = V (75V/) 



p 2 dp 



a 2 (x,p)p 



2 d£_ 
dp 



-(a 1 (x,p)-V)p 2 f + a s (x)p 4 f 



(3) 



Here the coefficients a\ and ai describe the second-order 
Fermi process and a s synchrotron losses. Equations (1) 
and (3) give the full advection-diffusion equation of cosmic 
particles in the diffusion approximation. 

The general system of SDEs describing the motion of 
a point X phase space can be written as: 



where Wx is a Wiener process: a stochastic diffusive pro- 
cess used to describe Brownian motion, with a conditional 
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probability which follows a Gaussian distribution. For ini- 
tial conditions given by Wx — Wx at t — to, we have at 
time t: 

{Wx) = W Xo , 
{(Wx-Wx ) 2 ) = t-t . (5) 

where ( ) means the average value. Ito ( |l95l| ) has shown 
that the distribution function describing the stochastic 
trajectories of the point X obeys the Fokker-Planck equa- 
tion if the coefficients A and B are identified with the dy- 
namic friction and diffusion coefficients of this equation. In 
the case of the Fokker-Planck equation (1), together with 
(2) or (3), the phase space is four-dimensional, X = (x,p), 
and expressions for the coefficients A and B are given 
by KA94. In this paper we shall restrict ourselves to the 
case of one space dimension, and include only the terms in 
Eq. (3) describing spatial diffusion and synchrotron losses, 
i.e., a± t 2 = 0, so that the phase space is two-dimensional: 
X = (x,p). Further simplifying to the case where the spa- 
tial diffusion coefficient is independent of both position 
and momentum and the loss rate a s is independent of po- 
sition, the set of SDE's (Q) reduces to 



d ln(p) 



da; = udt + V2k dW x , 
1 du 
3 da; 



a s p 



dt 



(G) 
(7) 



In this system, W x is a continuous but non-differentiable 
process, so that, as written, these equations do not exist in 
a strictly mathematical sense. To overcome this problem 
Ito ( 1951 ) (see also Gardiner 1983|) has defined stochastic 
differential integrals of the form J B(X(t), t)dWx- Ap- 
proximate solutions of the system of SDE's can be found 
by discretizing in time: 



Ax = x k+1 - x k 



\n(p k+ i/p k ) 



u(x k )(t k+ i - t k ) + V2k x 
[W(t k+1 ) - W(t k )} 
1 du 



a s Pk + r -r- 
3 dx 



(8) 

(tk+i-tk) (9) 



The term in brackets is the increment of the Wiener pro- 
cess, which is proportional to the square root of t k+ \ —t k = 
At: 



AW = £At 



1/2 



(10) 



where £ is a Gaussian distributed random number with 
zero mean and unit variance . This is the so-called Cauchy- 
Euler procedure (Gardiner 1983 ), as used by KA94. At 
each time step At, the change in x is made up of two 
parts: an advective (and deterministic) step 



ASadv = u(x k )At 

and a diffusive (stochastic) step 
Ax diS = £^2nAt . 



(11) 



(12) 



The approximate solutions converge to the exact solution 
of the SDE for At — ► 0. However, this scheme has the dis- 
advantage that At must be small enough to resolve the 
spatial structure in u(x). Denoting by x s hock the shortest 
lengthscale associated with u(x), KA94 found for a par- 
ticular example the requirement 



Axadv < Min (Axdiff, x 



shock J 



(13) 



Thus, although the diffusive step can be long compared to 
^shock, the advective step must resolve it, and the method 
is not appropriate for flow patterns containing sharp gra- 
dients (or discontinuities). 

2.2. Implicit Euler schemes 

Implicit methods, in which the increments arc expressed 
not just in terms of the solution at the start of a time 
step, but implicitly in terms of the solution at the end 
of it, are frequently effective in relieving time-step prob- 
lems, and h ave been discussed for SDE's by Smith & Gar- 
diner (|1989D . Their main advantage is that they yield sta- 
ble algorithms. However, the problem raised by the con- 
dition (13) is not one of instability, but accuracy. Further- 
more, it would appear that the advective, deterministic 
term is more sensitive to the problem than is the diffusive 
term. In view of this, we have chosen to test an algorithm 
in which, for the advective term, the coefficient is evalu- 
ated neither at the initial point of a time step (explicit) nor 
at the end point (fully implicit) but is integrated exactly 
over the step, using a linear interpolation of the trajec- 
tory. In this way, we may hope to account approximately 
for changes in the velocity u(x) which occur on very small 
length scales, unresolved by either Ax a dv, or Axdis . 

Replacing the advective terms in Eqs. (8) and (9) using 



Ax . 
At 



t 



(14) 



we find, neglecting for the moment the synchrotron losses, 



Ax 



dtu(x) + \/2K[W{t k -\ 



W(t k )} 



ti 

At 
Ax 



— / dxu(x) 



In 





+V2^ 


er) 


= 1 f 








At 



3Ax 



vv {lk+1) - 

+1 dt^ U 
dx 

[u (Xfc+l) 



W(t k )} 



(15) 



(16) 



which has the character of an implicit scheme, since the 
right-hand side of Eq. (15) is a function of x^+i . This tech- 
nique has been used i n app lications of SDE to other fields 
(see Kloden & Platen 1992| ) , but we are not aware of a de- 
tailed discussion of its properties. We show in the follow- 
ing that for the test cases we have examined, the scheme 
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yields accurate results when the condition (13) is replaced 
by the less restrictive one: 



(17) 



3. The test case of acceleration at a single shock 

We consider an infinite ID plasma in which particles prop- 
agate with diffusion coefficient kd- The flow velocity of the 
plasma is constant for x < — x s hock/2 (the upstream re- 
gion) and equal to u\ in the shock restframe. Similarly, 
in the downstream region, x > x s hock/2, the velocity is 
constant, U2 = u\/r, where r is the compression ratio of 
the shock. Position, time and momentum are normalised 
to the diffusion length and time scales upstream: (kd/ui, 
and Kv/ui respectively) and the injection momentum of 
a particle p;. Thus we define the following dimensionless 
variables and coefficients: 

x = xu 1 /k^ ) 

P = P/Pi 
i = tuf/nu 
u = u/ui 
a s = a s (x)/a sfi 

rx Bh ock/2 

t] = u\ I dx/nn 



J -a; S hock/2 
uia; s hock/KD 



(18) 



where a Si o = 3wf/(2KDp i ) controls the synchrotron losses, 
and r\ is the Peclet number, since the shock transition is 
confined to the region |a;| < 77/2. 

Using the scheme (15,16) we have simulated trajecto- 
ries in shocks of different r\. Each simulation runs with a 
given number of particles TV, injected at momentum p = 1, 
for a given computation time t m ax- Particles can escape 
at t < t max by crossing a boundary at |x| = L 3> 1. The 
particle distribution is measured at the shock front. Each 
crossing with an initial momentum p > 1 increases by 
unity the differential logarithmic distribution in the bin 
of momentum bracketing p. 

It is easily seen from Eq. (16) that for 77 -C 1 large 
values of Amp occur if the advective and diffusive steps 
almost cancel: Ai a( j v rs — Aidiff • In this case the result is 
particularly sensitive to our assumption that the trajec- 
tory between the points ik and Xk+i can be interpolated 
linearly. A partial solution to this problem is to reduce the 
time step which automatically decreases the probability of 
choosing a diffusive step which cancels the advective step. 
However, this procedure increases the computation time. 
The problem can be avoided completely by replacing the 
Wiener process by a different random process possessing 
the same mean and variance, but which effectively pre- 
vents the rare steps with Ai a d v ~ — A^diff- For a large 
number of trials, the random step produced by any such 



distribution tends to the Gaussian form of W. The sim- 
plest choice is to take £ = ±1 with equal probability for 
each sign. 

We first test this prescription in the most difficult case 
of infinitely thin shocks where rj — 0, in which case the 
profile shows a discontinuous change in velocity between 
the up and downstream regions. 

Figure [l] gives the stationary spectrum obtained with 
the modified scheme described above. We found over at 
least two decades of energy a stationary power-law distri- 
bution function of index s = 3r/(r — 1) = 4 in the case 
of strong shocks. Larger computation times are needed to 
reach the same accuracy at larger momentum. A \ 2 test 
against the analytic solution gives a value ~ 1. 




n 



i°g[p/p ] 



Fig. 1. Single shock stationary solution. The parameters 
are: compression ratio r = 4, maximum trajectory time 
t m ax = 3000, boundaries L = ±1000, time step At = 10" 3 . 
The x step is calculated implicitly using (15). Upper panel: 
comparison between numerical distribution and the analytic 
solution / tx p~ 4 (solid line). Lower panel: the distribution 
weighted with p 4 . 



We consider now the case of a shock with Peclet 
number r\ > 1. We compare our results with the semi- 
analytical derivation of the spectral index by Schneider 



& Kirk (1987). The simulations here use a linear velocity 
profile 



u(x) 



1 + r (r — l)x 



2r 



rjr 



(19) 



For larger Peclet number, the particles experience a 
smaller velocity jump at each step. The Fermi process is 
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thus less efficient, leading to softer stationary distribu- 
tions. For a Peclet number of 77 = 10 we get an index of 
~ 5.8 (see Fig. 0) in good agreements with the results of 



Schneider & Kirk (1987). The discrepancies seen at small 
momenta are not statistical errors, but arise because our 
method assumes zero particle flux in space at the injection 
point, whereas the scale- independent power-law distribu- 
tion implies a finite diffusive flux. This 'transient' effect 
disappears at momenta slightly above that of injection. 



where the first-order Fermi term (Fi) and the loss term 
L s are 



F\ = exp 



and 



Au 
3A5 



■At 



At 

'A~x Pk 



exp 



Ai 
3Ar" 



At dx' . 



(22) 



(23) 




This solution is exact for the linearly interpolated tra- 
jectory (14). Having determined Xk+i from Eq. (15), the 
new value of the momentum is given by Eq. (21) with 

X = Xk+l- 

In Fig. ^, we show the results of this implicit scheme. 
Following Webb et al. ( |1984| ) and KA94, we define up- 
stream and downstream synchrotron coefficients (a s i iS 2), 
given by a S i = u 2 / (4ka s pi). This gives for the charac- 
teristic momentum = 4/(g/a si + (q — 3)/a S 2), with 
q = 3r/(r — 1). For a compression ratio r = 4, we have 
a s i — 10 and a S 2 — 1, and the distribution cuts-off for 
p > p* — 2.9 (i.e., log(p) > 0.46). The solution showed in 
figure 3 is in good agreement with the analytical result of 
Webb et al. 



(1984) 



1.5 



log P/Pol 



Fig. 2. Upper panel: the distribution for same parameters as 
in Fig. (Q), but with the Peclet number f] = 10 and time step 
At = 10~ 3 , compared to the analytical solution / oc p -5,8 . 
Lower panel: the spectrum for single shock weighted by p 5 ' 8 . 



3.1. Synchrotron losses 

Synchrotron losses in the diffusive shock acceleration 
process have been investigated analytically by Webb et 
al. (1984). Their main effect is to soften the spectrum at 
momenta greater than p* where the loss rate equals the ac- 
celeration rate. The inclusion of loss terms in the scheme 
modifies the way the momentum gain is calculated. Re- 
turning to Eq. (7), and using the linear interpolation of 
the trajectory (14) we arrive at the ordinary differential 
equation 



-V- 



(20) 



dp 1 d " Ai 

dx sP Ax 3 dx Ax 

For the initial condition p = p^ at x — ik (i.e., t — t^) the 
solution is 






log[p/Po! 




Fig. 3. Upper panel: the distribution for acceleration at a single 
shock front including synchrotron losses for a s i = 10, a S 2 = 1 
and r = 4. The thin solid lin e is ad apted from the analytical 
solution given by Webb et al. ( 1984 ) . The maximum time for a 
single trajectory is t ma x = 3000, and At = 10~ 3 . Lower panel: 

i, the solid line 



In (p/pk) = - In (Fi + L s ) 



(21) 



, , , lx = 3000, and At = 10 
the distribution weighted by p 4 . In each cas 
corresponds to the loss free solution. 



(i 
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The implicit scheme is, of course, much faster than the 
explicit scheme for flows with small Peclet numbers, since 
the latter are accurate only when At <C rj, whereas for the 
former At <C 1 suffices (see Table Q). However, it is also 
interesting to compare the implicit and explicit schemes 
in flows with large Peclet numbers, where both should be 
accurate. In this case, the explicit and implicit schemes 
give similar results for single shocks for At ~ 10~ 3 — 10~ 2 . 
For larger time steps the particles tend to be advected 
prematurely from the acceleration zone and both schemes 
show an unphysical softening of the spectrum. For a given 
At, the explicit scheme is faster by a factor ~ 5/3, since 
all quantities entering into the computation of the position 
and momentum increments are calculated only once per 
time step. 



Table 1. Summary of the results of the explicit and implicit 
schemes. The power-law index of the distribution is compared 
with the v alue given by an analytic calculation (Schneider & 



Kirk 
of ~ 



1987 ). The typical error on the index estimates is of order 



few%. 



V 


Analytic 


At 


implicit 


explicit 


IO" 12 


s = 4.0 


10~ 3 


s=4.0 




210" 3 


s = 4.0 


10" 3 


s=4.01 


s=3.98 


21(T 3 


s = 4.0 


10~ 2 


s=3.99 


s=4.55 


21(T 3 


s = 4.0 


ht 1 


s=4.04 




1 


s = 4.18 


10~ 3 


s=418 


s=4.19 


1 


s = 4.18 


10~ 2 


s=4.18 


s=4.19 


1 


s = 4.18 


io- 1 


b=4.21 


s=4.21 


10 


s = 5.80 


IO" 3 


s=5.80 


s=5.79 


10 


s = 5.80 


10~ 2 


s=5.81 


s=5.82 


10 


s = 5.80 


IO" 1 


s=5.85 


s=5.92 



4. Acceleration at multiple shocks 

The subject of multiple shock acceleration has been exten- 
sively investigated both analytically and numerically over 
the past few years (Spruit 1988, Achterberg 1990, Schnei- 



der |1993| , Pope & Melrose 



1993 



Melrose & Crouch 1997) 



Analytic solutions of the diffusion-advection equation in- 
cluding synchrotron losses can be derived if we assume 
that the time spent by a fluid clement between two con- 
secutive shocks is much longer than the acceleration time 



at a single shock problem (see Schneider 1993). Another 
way of formulating this condition is that the first-order 
Fermi acceleration process should be much faster than all 
other processes, such as escape, decompression and losses. 
With this hypothesis, a final power-law index can be cal- 
culated which takes account of the different generations of 
particles accelerated at individual shocks. At sufficiently 
high momentum, this approximation fails, since the time- 
scale of the losses must eventually become comparable to 
the acceleration time-scale. The effect of multiple shocks 
is to increase the acceleration efficiency, by reducing the 



effective escape rate. In a purely ID system, the escape 
rate is formally zero, and the distribution functions tends 
to f (x p~ 3 . As shown, using spatially averaged equa- 
tions, by Kardashev ( 1962| ) and by Schlickeiser (1984), in 
the presence of losses, this spectrum extends up to mo- 
mentum values where the loss effect becomes dominant, 
and the spectrum piles up at a momentum where the ef- 
fective acceleration time equ als th e loss time. More re- 
centl y, Pr otheroe & Stanev ( 1999 ) (see, however, Drury 
et al 1999 ) have proposed an alternative method of com- 
puting the cut-off and pile-up effects in the high energy 
particle spectrum with various energy dependent diffusion 
coefficients. They present both an analytical model and a 
conventional Monte-Carlo simulation (of the same kind as 
described in the introduction) and show that the effect of 
the Klein-Nishina regime of inverse Compton losses may 
modify the results obtained with continuous (synchrotron 
and inverse Compton in the Thomson regime) losses. In 
this paper we include neither non-continuous energy losses 
nor energy dependent diffusion, but postpone an investi- 
gation of these effects to future work. 

In general for multiple shock system, the picture is 
somewhat more complicated than depicted in previous an- 
alytical or semi-analytical models, since the shocks may 
be so close together that the acceleration time at a sin- 
gle shock is not short compared to the flow time between 
the shocks. Also, at high momenta, the synchrotron loss 
time-scale must be compared not only to the acceleration 
time at a single shock, but also to the flow time between 
the shocks. 

To investigate these situations, we consider a simple 
periodic pattern (of period L) including shock fronts and 
re-expansion regions (see Fig. [|). The flow speed of the 
pattern [x € [-L/2, L/2]) can then be written as 



u(x) = u\ + (ui 



, 2 j ^ +Xs/2 ,xG [-L/2,-x a /2] 



L i s 



~r\ ~ /- _n5 + x s /2 

U{X) = U\ — (Ui — U2) , X <E 



~x s /2, Xs/2\ , 



u{x) = u 2 + (ui - u 2 ) X Xs ^ , x £ [x s /2, L/2] (24) 
L x s 

where 112 = iii/r is the flow speed on the downstream side 
of the shock. 

There are two free parameters (for a given compres- 
sion ratio) in this problem: the advection time t a dv, which 
is the (dimensionless) time for the fluid to flow through 
one wavelength of the pattern and is numerically equal to 
the dimensionless length L, and the loss strength given by 
a s , which effectively defines a momentum scale, since the 
(dimensionless) synchrotron loss time at momentum p is 



3a s p 



(25) 



To these we add a third: an escape time t osc , assumed in- 
dependent of x and p. This can be understood as a crude 
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Fig. 4. Geometry of the multiple shock system - a single period 
of the pattern is drawn in bold. 



way of incorporating 2 or 3-dimensional effects into our 1- 
dimensional simulation, since only in 1 dimension are the 
particles unable to leave the train of shock fronts. The in- 
clusion of escape effects in the SDE system can be effected 
by rescaling the number of particles that have crossed the 
shock with at time t by the factor exp(— t/t eBC ). 

There are then two important characteristic times 
which are important for a discussion of the solutions: 

— the single shock acceleration time ts , which is the char- 
acteristic time of momentum increase at an isolated 
shock front: 

t s = 3ri±; , (26) 



1 



(Lagage & Cesarsky 1981) with k = k~i = k~2 = 1 



— the effective multiple shock acceleration time f]yt which 
is controlled by both the advective time and the sin- 
gle shock acceleration time and is inferred from our 
numerical results. 

We present, in the following, results for a shock system 
described by Eq. (24) with a Peclet number rj = 0. We 
adopt as typical values of the parameters: compression 
ratio r = 4, escape time t esc = 10 3 , an inter-shock distance 
L = 20, and, in the case of synchrotron losses, a s = 6.5 x 
10~ 6 (chosen to give a peak momentum of 10 3 ). All the 
simulations are run for a time of i max = 5t esc and with a 
time step At = 10~ 2 . 

4-1- Multiple shock effect: the case without losses 

We first consider the case where losses are inefficient 
(5 S — > 0). The results are shown in Fig. [|a. The sta- 
tionary solution at each shock front fix = nL,p) (with 
n an integer) is close to p~ 3 . This confirms the multi- 
ple shock effect as an efficient way of producing higher 
energy particles and harder spectra than isolated shocks. 



In fact, the stationary index is not exactly 3 but slightly 
steeper, owing to the finite escape time. Using the general 
relation between the power-law index of accelerated parti- 
cles and the escape and acceleration times: / oc p~ s , with 
s = 3 + f acc /fesc (Kirk et al. 1994) we find for the effective 



acceleration time in this particular multiple shock system 
<m ~ 10 2 . Thus, t M is (= L) in this case. 

For momenta lower than that of injection, p < 1, 
t he m ethod of Schneider ( 1993 ), and Melrose & Crouch 
( 1997 ), gives a continuation of the power-law down to 
p ~ 1/r. This is clearly an artifact of the assumed separa- 
tion of the acceleration and expansion processes. It does 
not appear in our simulations, which show instead a hard- 
ening to lower frequencies starting at the point p = 1. 
However, because of the transients associated with our 
method close to the injection momentum, this effect may 
also be an artifact. 

At still lower momentum values, the analytical solution 
f(p) oc p x for an escape probability independent of mo- 



mentum is given by (see equation 4.5 in Schneider 1993) 



(3r/(r-l)-3)(r 1+A / 3 ) 
3r/(r- 1) + A 



= 1 



^adv 



(27) 



The stationary index is A = for t csc 3> £ a dv, but for 
lower values of the ratio of escape to multiple shock ac- 
celeration time, the spectrum hardens, and typically for a 
ratio tadvAcsc ~ 0.1, we get A ~ 0.2, in good agreement 
(given the accuracy of the time-scales) with the simula- 
tions; for p < 1 f(p) ~ p 3 down to p = 0.1. 

4.1.1. Variations of the escape time 

If the inter-shock distance L is kept constant, increasing 
(decreasing) the escape time leads to a harder (softer) sta- 
tionary spectra. This is clearly seen in Fig. ^b where we 
have reduced the escape time by a factor of 2 compared 
to the fiducial case of Fig. |[ The resulting spectrum has 
an index of 3.2, which again gives tyi — fewxiadv For mo- 
mentum p < 1, from Eq. (p7|), the relation t csc ~ 5?m gives 
a spectrum with A ~ 0.4. We obtained down to p = 0.1 
an index of ~ 0.5. If t csc < t&dv, the particles escape the 
system before being advected to the next shock. The sta- 
tionary solution tends to the single shock result of Sect. 2.1 
without losses, which is a power-law spectrum with an in- 
dex of 4 (for r = 4). 

4.1.2. Variations of the inter-shock distance 

We now keep t csc constant and equal to the fiducial value 
of 100. For intershock distances L = < ac i v S> 100, no mul- 
tiple shock effect is seen (Fig. ||c) . As in the case of short 
escape time (Fig. ||b), the stationary solution tends to that 
from a single shock. As i a d v decreases, the spectrum hard- 
ens and the spectral index can take all values between 3 
and 4 (for r = 4). This is consistent with our finding that 
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tyi few x i a dv- If the inter-shock distance is reduced 
to below one diffusion length L < 1, the assumptions nor- 
mally used to derive the basic transport equation (Skilling 



of this hump on the low momentum side is determined 



by the momentum at which t c 



t S yn, where escape in- 



1975) are violated. We have not investigated this regime; 



although it potentially interesting, since if the distribution 
remains almost isotropic, the diffusion approximation may 
in fact remain adequate. For low momentum p < 1, the 
spectrum steepens with decreasing L and tends to p°. 





2 4 

l°g[p/p„] 



2 4 

log[p/p„] 



Fig. 5. Differential logarithm spectrum produced by an en- 
semble of identical shocks. Upper left panel (a): our refering 
case with tM/tesc ~ 0.1 (see text for parameters). Lower left 
panel (b): the case with different escape times, in solid line 
iosc = 1000, in dotted line t csc = 500, and in short-dashed 
line t csc = 100. Lower right panel (c): the case of different in- 
ter-shock distances, from the upper to the lower curve we have 
dashed-dotted line L — 10, solid line L = 20 our fiducial case, 
long-dashed line L = 50, dotted line L — 100, short-dashed 
line L = 10 4 . The power-law spectrum corresponds to the sin- 
gle shock solution. 



4-2. Multiple shock effect: the case with synchrotron losses 

The stationary spectrum computed with the fiducial pa- 
rameters is given in Fig. |^a. The inclusion of losses cre- 
ates a pile-up at a momentum where losses equal gains 
i.e., <m = t syn . From the simulation, we find this occurs at 
roughly p — 10 3 , which implies tyi ~ 100 close to 5 x t a d v - 
This provides an independent check on the estimate made 
from the results presented in Figs. |^b and ||c. The width 



tervenes to overwhelm the effects of synchrotron losses. 
Since t esc = 50t a d v ~ IO^m, the hump makes its appear- 
ance about one order of magnitude before it peaks. 

At momenta higher than the peak of the hump, the 
losses dominate over all other processes, since there t syn < 
tadv < £m, and, in this example, the acceleration at a sin- 
gle shock is approximately equal to the advection time 
(see Eq. 26). Consequently, the spectrum cuts off expo- 
nentially. This conclusion is consistent with the results of 



Melrose & Crouch (1997) found using an iterative method 



in which the synchrotron cut-off is estimated in advance. 

4.2.1. Variations of the loss rate 

A decrease (increase) of d s with L and t osc kept unchanged, 
leads to stronger (weaker) losses, and lower (higher) mo- 
menta at the peak of the pile-up. This is confirmed in 
Fig. ^|b, which also demonstrates that the peak momen- 
tum is directly proportional to the loss time t syn . The 
width of the hump is unaffected, since its lower bound also 
moves in proportion to t syn . Radiative losses are unimpor- 
tant below the injection momentum, except when the loss 
rate is so strong (a s ~ 1) that particle are prevented from 
being accelerated at all. 

4.2.2. Effects of the other parameters 

We stress here the effects of both escape and adiabatic 
losses on the shape of the pile-up. 

— Keeping the inter-shock distance L = i a d v constant 
and allowing t csc to vary causes a change in both 
the spectral slope at low momenta, where losses are 
unimportant, and a change in the width of the pile- 
up Fig. ^c. As t csc increases, the low momentum spec- 
trum hardens, and the width of the hump increases, 
as it shortens, the spectrum softens, and the hump be- 
comes less prominent, until it disappears completely 
once £ C sc < iadv At this point, the power-law index 
of the low momentum spectrum is approximately 4, in 
agreement with the findings of Kardashev ( 1962| ) . 

— Keeping the escape time constant, but varying the ad- 
vection time (and inter-shock distance) enables one to 
distinguish the regimes of single and multiple shock ac- 
celeration. For larger values of t a d v , than our fiducial 
case, four regions in momentum space can be found. 
Starting at low momentum, these are (Fig. ^|d) 

1. For f syn > t csc , acceleration proceeds by the multi- 
ple shock process without losses. 

2. For t esc > t syn > t a dv the spectrum is affected by 
losses, and starts to form a pile-up. Acceleration at 
multiple shocks takes place. 

3. For t a dv > ^syn > ts particles are prevented from 
reaching the next shock before cooling. An interest- 
ing phenomenon appears in this case: the spectrum 
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shows a power-law index appropriate to accelera- 
tion at a single shock front (see the dotted line) 
4. For is > tsyn losses dominate over all other pro- 
cesses and the spectrum cuts off. 
In Fig. ^|d all four regions can be distinguished in the 
case L = 100. At higher L, escape prevents multiple 
shock acceleration, and at lower values the shocks are 
too close together to allow the emergence of the single- 
shock power-law spectrum. These effects do not appear 
in previous works, which had to assume a separation in 
the time-scales of the different. The high momentum 
power-law tail due to the single shock process may be 
important and contribute a substantial fraction of the 
total pattern luminosity. 



4-3. Synchrotron spectrum 

The optically thin synchrotron emission produced by the 
particle distribution within a shock pattern described by 
Fig. ||, assuming constant magnetic field (i.e., parallel 
shocks) is given by 



S v oc 



L/2 



L/2 



P{v,p)f{p,x)p dpdx 



(28) 



where, for a relativistic particle of charge q, mass m, and 
pitch-angle 



P[y,p) 



q 3 B sintf 



2tt 



(29) 



c J v I i/ c 




2 4 

l°g[p/PoJ 




2 

l°g[p/Pol 



The characteristic frequency depends on p as v c {p) = 
3u> c sin t? p 2 /(2m 2 c 2 ), where the cyclotron frequency is 
lu c = qB /mc. The integration over the particle position is 
done by computing the number of particles at a given mo- 
mentum p at each time step inside the pattern. This gives 

L 12 

the quantity N(p) — f_ L / 2 f(P->%)P 2 dp contributing at a 
given frequency by the factor P{v,p) to the synchrotron 
spectrum produced by the pattern. Figure [?] shows the re- 
sulting synchrotron spectrum. 



Fig. 6. Differential logarithm spectrum for multiple shock ac- 
celeration with synchrotron losses. Upper left panel (a) our fidu- 
cial case, see text for the parameters. Upper right panel (b): the 
case with different synchrotron loss rates, the solid line corre- 
sponds to our refering case, the short dashed and dotted lines 
have a loss rate divided by 10 and 100, and the long dashed has 
a loss rate multiplied by 10. Lower left panel (c): the case with 
different escape time, from the upper to lower curve: t csc = 
210 3 , 10 3 our refering case, 500, and 100. Lower right panel 
(d): the case with different inter-shock distances, the distri- 
bution correponding to our refering value L = 20 is in solid 
line. In dotted-dashed line L = 10, in long-dashed line L = 50, 
in dotted line L = 100, in short-dashed line L = 10 4 . The 
power-law spectrum corresponds to the single shock solution 
with ts/iosc ~ 210" 2 . 




iog|>] 



Fig. 7. Synchrotron spectrum produced by a periodic pattern. 
The solid lines line shows the synchrotron radiations produced 
if the loss rate is one hundred times stronger (see the solution in 
dotted line of figure 6 b). We compare also the low momentum 
slope with a 1/3 power- law spectrum is dashed line. The dotted 
line curve is produced for an escape time a factor two larger 
than in our refering case (see the solution in long-dashed curve 
of figure 6 c). 
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The results, shown in Fig. [7], display a great variety 
of synchrotron spectra. For large ratios i e sc/iadv the spec- 
trum is flat at low frequency, and extends over a range 
which depends on the relative strengths of loss and accel- 
eration, which can be quite sufficient to explain the flat 
radio spectra observed in radio-loud quasars. Such an ex- 
planation provides an alternative to the inhomogeneous 
self-absorbed models usually advanced as the origin of fiat 
radio Quasar spectra. Figure (?] also shows that inverted 
spectra over 2-3 decades are also quite possible, where the 
cooling time is still shorter than the escape time. In fact 
the pile-up effect, as stressed by Melrose ( |l996| ), is essen- 
tial to explain the inverted spectra observed in galactic 
centre sources such as Sagittarus A* (Beckert ct al. 1996 ), 
and also may be associated with the flaring states of some 
extragalactic sources, such as the recent X-ray bursts of 
Mrk 501 (Pian et al. 1998). Such flares can arise from a 



variation of either the inter-shock distance or the escape 
time from the system. A change in the magnetic field 
stength does not mimick the multi-wavelength behaviour 



adequately (Mastichiadis & Kirk 1997) 



5. Conclusions and outlook 



Achterberg & Kriills (1992) and Kriills & Achter 



berg (1994, KA94) have demonstrated the usefulness of 



the stochastic differential equation approach as an effi- 
cient tool for investigating particle acceleration at shock 
fronts. They found approximate solutions of the diffusion- 
advection equations in several different physical situa- 
tions including second order Fermi acceleration and mo- 
mentum dependent spatial diffusion coefficients. However, 
their computational approach is limited to spatially re- 
solved shock structures, and needs unrealistically short 
time steps for thin shocks. We have presented an improved 
scheme in which the particle momentum gain during a 
time step depends on both the initial and final positions 
of the particle. The scheme reproduces analytical results 
derived for shock thicknesses much lower than a typical 
diffusive length, without imposing an excessive require- 
ment on the time step. We have applied our procedure 
to a system of multiple identical shocks. The p~ 3 signa- 
ture is obtained when the escape time is much larger than 
the multiple shock acceleration time. Inclusion of losses 
leads to a pile-up effect where the acceleration and loss 
rates are equal, as has been found analytically. The pile- 
up is present only for a restricted range of parameters: 
when the escape time exceeds the time for the plasma to 
move from one shock to the next, but is short enough not 
to permit synchrotron cooling of the lowest energy parti- 
cles. We have shown that parameter ranges can be found 
in which the spectrum simultaneously displays the power 
law index characteristic of multiple shock acceleration (at 
low momenta) and the index appropriate to single shock 
acceleration at high momenta. 



The synchrotron spectrum produced in each flow pat- 
tern can show flat or inverted spectra depending on the 
inter-shock distance, the loss strength and the ratio of 
escape to acceleration time. Multiple shock acceleration 
may explain the radio spectrum from flat radio Quasars, 
the radio to IR spectrum of galactic sources such as Sagit- 
tarus A*, and even hard X-ray spectra observed in high 
states of the BL Lac Markarian 501. 

Several extensions of this work are possible, such as the 
consideration of non-periodic flows and non-stationary in- 
jection. It also seems feasible to extend the ID simulations 
to more complex multi-dimensional flows, such as numer- 
ical simulations of jets. 
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